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I .   INTRODUCTION 

A.   BACKGROUND 

In  analyzing  ocean  acoustic  pulse-propagation  problems, 
accurate  phase  calculations  must  be  performed  in  order  to 
predict  the  acoustic  signal  at  the  receiver.  Linear  systems 
theory  provides  analytical  expressions  for  analyzing  the 
propagating  acoustic  field  [Refs.  1-4].  These  well-known 
equations  form  the  basis  for  the  FORTRAN  programs  developed  in 
this  thesis. 

Program  input  is  depth  versus  sound-speed  data  pairs  where 
the  speed  of  sound  is  an  arbitrary  function  of  depth.  The 
ocean  is  viewed  as  a  linear,  time-variant,  space-variant 
filter.  The  WKB  approximation  can  specify  this  filter's 
transfer  function  [Refs.  1-4].  For  the  arbitrary  sound- 
speed  profile,  no  exact  transfer  function  exists.  The 
transfer  function  requires  a  method  for  evaluating  signal 
phase  at  the  receiver. 

Two  phase  evaluation  methods  are  presented  and  contrasted. 
Each  method  calculates  the  phase  of  the  acoustic  signal  for  a 
specified  horizontal  range.  The  first  method  overlays  a 
piecewise  linear  profile  on  the  input  sample  values  for  the 
speed  of  sound  at  various  depths.  The  theory  of  ray  acoustics 
is  used  to  calculate  travel  time  and  phase.  The  second  method 
fits  a  smooth  cubic  spline  curve  to  the  input  samples.   This 


method  solves  the  propagation  problem  using  a  system  of  three, 
first-order  differential  equations  [Ref.  5].  Both  methods 
make  phase  calculations  using  signal  travel  times. 
Additionally,  each  method  is  capable  of  identifying  eigenrays 
or  rays  that  directly  connect  the  signal  source  to  the 
receiver. 

Another  method  is  introduced  to  validate  results.  Direct 
numerical  integration  is  performed  to  calculate  phase  for  a 
single  gradient,  linear,  sound-speed  profile.  A  separate 
FORTRAN  program  implements  this  numerical  integration  routine 
providing  a  totally  independent  verification. 

An  overview  of  the  theory  behind  each  method  is  presented 
in  Chapter  II.  Computer  simulation  results  are  presented  in 
Chapter  III  for  various  input  sound-speed  profiles.  Tabular 
results  of  phase  calculations  are  presented  for  each  of  the 
three  analysis  methods.  Using  input  from  the  piecewise  linear 
and  the  cubic  spline/differential  equation  solution  methods, 
ray  traces  of  the  propagating  field  are  presented  as  a  visual 
aid  in  interpreting  the  results. 


II.   THEORETICAL  BACKGROUND  FOR  THE 
EVALUATION  OF  PHASE  INTEGRALS 


A.   UNDERWATER  ACOUSTIC  PULSE  PROPAGATION 

Ocean  acoustic  pulse-propagation  models  can  be  derived  by 
using  the  principles  of  linear,  time-variant,  space-variant, 
systems  theory  and  the  physics  of  wave  propagation  in 
inhomogeneous  media.  Linear  systems  theory  allows  for  the 
development  of  an  ocean  medium  transfer  function.  An  ocean 
medium  transfer  function  that  is  based  on  the  WKB 
approximation  has  been  derived  and  is  given  by  References  1 
and  2  as  follows: 


ff<f,fP,y„;y>  =        A        e~lJ>,dc  e<*".<™>         (2.i) 
/lMy)T 


where 


f     is  the  frequency  in  Hz 

fr    is  the  radial,  spatial  frequency  in 
cycles  per  meter 

y0    is  the  source  depth  in  meters 


A  =  j  y/1  My0)  I  (2.2) 


47Tfy 


k^(y)  is  the  propagation  vector  component  in 
the  Y  direction  with  units  of  radians 
per  meter  and  is  given  by 


My>  = 


±    27r{[f/c(y)]2   -    fr2)h    ,fr   <    f/c(y)  (2.3) 

[T   j27T{fr2   -    [f/c(y)]2}'/2    ,fr   >    f/c(y)  (2.4) 


c(y)  is  the  sound-speed  expressed  as  a  function 
of  depth  with  units  of  meters  per  second 

f   is  the  transmitted  (input)  spatial  frequency 

in  the  Y  direction  at  the  source  in  cycles  per 
meter  and  is  given  by 


fy    = 


±    [(f/c0)2    -    fr2]*    ,    fr    <    f/c0  (2.5) 

=F       [fr2    -    (f/c0)Y2    ,     fr    >    f/c0  (2.6) 


and 


c0   is  the  speed  of  sound  in  meters  per  second  at 
the  source  depth  y0,  that  is,  c0  =  c(y0); 

The  plus  (minus)  sign  in  Eqs.  (2.3)  and  (2.5)  is  chosen 
whenever  y  -  y0  >  0  (y  -  y0  <  0)  .  The  minus  (plus)  sign  in 
Eqs.  (2.4)  and  (2.6)  corresponds  to  the  plus  (minus)  sign  in 
Eqs.  (2.3)  and  (2.5).  This  thesis  deals  only  with  propagating 
waves.  Therefore,  Eqs.  (2.4)  and  (2.6)  representing  the 
generation  of  evanescent  waves  (i.e.,  decaying  exponentials), 
will  not  be  used. 

The  evaluation  of  the  phase  integral  in  Eq.  (2.1) ,  namely, 


6y(y)  =fyk  (C)dC  (2.7) 

y     Jyo  y 


has  been  attempted  in  past  studies  in  which  solution 
techniques  included  direct  integration  and  binomial  expansions 
[Ref.  3]. 

Computer  programs  using  these  techniques  proved  to  be  very 
expensive  in  terms  of  CPU  time.  This  thesis  compares  four 
different  methods  for  evaluating  the  phase  integral.  The  four 
methods  are  presented  in  the  following  sections. 

B.   DIRECT  INTEGRATION 

Computer  code  was  developed  to  directly  evaluate  the  phase 
integral  given  by  (Eq.  2.7)  using  numerical  integration 
routines  from  the  International  Math  and  Statistics  Library 
(IMSL) .  Direct  integration  is  a  time-consuming  technique  that 
was  applied  only  under  the  following  constraints: 

-  free-space  acoustic  propagation 

-  linear  sound-speed  profile  with  a  single  gradient  g 
in  inverse  seconds  and 

-  propagating  signals  that  have  not  passed  through  a 
turning  point  in  the  medium. 

Expressions  for  the  propagation  vector  component  k^(y)  of 

the  phase  integral  are  given  by  Eqs.  (2.3)  and  (2.4).   For  a 

linear  sound-speed  profile  with  a  single  gradient,   the 

dependence  of  c(y)  on  depth  y  is  given  by 

C(Y)   =  C(y0)  +  g(y  -  y0)  (2.8) 


where 


c(y0)  is  the  sound-speed  at  the  source  depth 
y0  meters 


g     is  the  sound-speed  gradient,  and 

y     is  the  desired  depth  for  sound-speed. 

Since  the  direct  integration  approach  is  limited  to  free- 
space,  single  gradient,  propagation  problems,  it  is  used  only 
to  validate  the  results  of  other  solution  techniques  for 
several  simple  test  cases.  The  results  will  be  compared 
against  the  next  two  methods  to  be  presented.  These  methods 
will  apply  the  theory  of  ray  acoustics  [Ref.  4  &  5]  to  find 
travel  time.  Phase  calculations  will  be  based  on  travel  time 
calculations . 

C.   PHASE  CALCULATIONS  BASED  ON  PIECEWISE  LINEAR  SOUND-SPEED 
PROFILES 

A  FORTRAN  program  developed  by  Lim  [Ref.  5]  was  the  ray 
acoustics  propagation  code  that  was  used  to  draw  ray  trace 
plots  and  to  calculate  travel  time,  in  addition  to  the  depth, 
angle  of  propagation  and  path  length  along  a  ray  path  as  a 
function  of  horizontal  range.  This  method  applies  ray 
acoustics  to  a  piecewise  linear  sound-speed  model  of  the  ocean 
medium.  Based  on  the  input  sound-speed  versus  depth  data 
pairs,  layers  are  defined  in  the  ocean  medium.  Each  layer  has 
an  upper  and  lower  boundary  at  specific  depths.  A  constant 
sound-speed  gradient  g  is  calculated  for  each  layer.  The 
sound-speed  for  any  desired  depth  is  computed  using  Eq.  (2.8) 
with  the  appropriate  gradient. 

In  a  constant  gradient  medium,  ray  acoustics  theory  allows 
calculation  of  travel  time  with  closed-form  equations.  For  an 
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incremental  increase  in  horizontal  range  from  the  source,  the 
angle  of  arrival  for  the  propagating  ray  can  be  shown  to  be 
[Ref.  5]. 


-1 


/3(y)      =  cos       [cos/30   -  b*g(y)    *   rngstp] 


(2.9) 


where 


(3Q  is  the  ray  launch  angle, 

g(y)   is  the  gradient  in  the  layer  at  depth  y,  and 

rngstp  is  the  incremental  increase  in  horizontal 
range . 

The  ray  parameter  b  in  Eg.  (2.9)  is  given  by  [Refs.  4  &  5]. 

b  =  sin/30  /  c(y0)  (2.10) 

Knowing  the  arrival  angle  )3(y),  the  ray  depth  y  is  given  by 
[Refs.  4  &  5] 


Y  =  Yo  +  — 


sinff(y) 
sin/?(yn) 


-1 


(2.11) 


The  arrival  angle  also  allows  calculation  of  travel  time  from 
[Refs.  4  &  5] 


t   =  —In 

g 


tan[/3(y)/2] 
tan[/3(y0)/2 


(2.12) 


Once  the  travel  time  (in  seconds)  is  found  for  the  desired 
horizontal  range,  total  phase  is  given  by 


9T  =  27TfT 


(2.13) 


The  phase  integral  given  by  Eg.  (2.7)  represents  a  phase 
change  in  the  depth  (y)  direction.  The  constant  value  of  the 


propagation  vector  component  in  the  radial  direction  allows 
for  easy  calculation  of  phase  change  in  the  radial  direction. 
This  value  is  given  by 

9R  =  27rfr  *  hrzrng  (2.14) 

where  hrzrng  is  the  total  horizontal  range  traveled. 
The  phase  in  the  depth  direction  is 

ey  "  eT  "  eR-  (2.15) 

Phase  expressed  in  radians  is  a  modulo  2n  function.  The 
proper  solution  for  the  phase  integral  Eq.  (2.15)  must  be 
expressed  as  a  modulo  2ir    function: 

6y  =  modulo  (9y  ,  27T  )  .  (2.16) 

This  method  (referred  to  as  method  1  in  Ref .  5)  is  capable 
of  analyzing  propagating  rays  at  any  horizontal  range  in 
either  free  space  or  a  bounded  medium. 

D.   PHASE  CALCULATIONS  BASED  ON  AKIMA  CUBIC  SPLINES  AND 
ORDINARY  DIFFERENTIAL  EQUATIONS 

This  application  of  ray  acoustics  applies  a  more 

sophisticated  treatment  to  the  depth  versus  sound-speed  data 

pairs  sampled  from  the  ocean  medium.   The  data  pairs  are  used 

to  form  a  smooth  Akima  cubic  spline  sound-speed  profile. 

Splines  offer  the  advantage  of  using  all  data  points  in 

generating  a  profile,  and  they  place  no  restrictions  on 

spacing  between  data  points.   The  Akima  version  of  the  cubic 

spline  was  chosen  for  its  excellent  ability  to  combat  wiggles 

in  the  profile,  that  is,  it  suppresses  oscillations  that  would 


cause  overshoots  and  undershoots  in  the  sound-speed  versus 
depth  profile. 

This  method  uses  ray  acoustics  theory  to  generate  a  system 
of  three,  first-order  differential  equations  [Ref.  5]: 

y,  =  y2  (2.17) 

y2  =  -n (2.18) 

b2c\y,) 


and 


Y3  =  z (2.19) 

bc2(y,) 


where 

c(y.,)  is  the  derivative  of  the  sound-speed  with 
respect  to  depth  at  depth  y1 , 

y1     is  the  ray  depth, 

y2    is  the  cotangent  of  the  ray's  arrival  angle 
0(y)  ,  and 

y3    is  the  travel  time  of  the  ray. 

Once  solved,  this  system  of  differential  equations  allows 

phase   to   be   calculated   from  travel   time.     The   phase 

calculation  is  performed  exactly  as  shown  for  the  piecewise 

linear  sound-speed  profile. 

Like    the    previous    method,    the    Akima    cubic 

spline/differential  equation  solution  can  be  used  on  free 

space  and  bounded  media  problems. 
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E.   EIGENRAYS 

Eigenrays  are  propagating  rays  that  exactly  connect  the 
sound  source  to  the  receiver.  The  FORTRAN  propagation  code 
developed  can  search  for  and  identify  eigenrays.  The  input 
required  is 

-  depth  versus  sound-speed  data  pairs, 

-  source  depth  in  meters, 

-  receiver  depth  and  range  in  meters, 

-  angle  step  in  degrees  between  possible  eigenrays  to 
be  evaluated,  and 

-  allowed  depth  error  denoted  yerror  in  meters. 

Rays  passing  within  the  allowed  error  or  tolerance  yerror  of  the 
receiver  are  identified  as  eigenrays.   The  eigenray  mode  can 
employ  either  the  piecewise  linear  sound-speed  profile  or  the 
Akima  cubic  spline/differential  equation  method  for  ray 
propagation. 

The  eigenrays  are  found  by  trial-and-error.  This  search 
method  was  chosen  after  the  IMSL  DBVPMS  program  failed  to 
solve  the  problem.  The  DBVPMS  program  is  a  differential 
equation  solver  that  was  applied  to  the  system  of  differential 
equations  given  by  Eqs.  (2.17)  through  (2.19).  The  routine 
uses  the  shooting  method  to  find  eigenray  solutions  to 
boundary  value  problems.  It  was  unable  to  converge  to  a 
solution.  Solutions  for  the  acoustic  problems  investigated 
are  difficult  for  this  algorithm  because  of 
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-  long  propagation  distances, 

-  the   inhomogeneous   ocean  presents  a   continuously 
varying  medium  with  discrete  boundaries,  and 

-  the  system  of  differential  equations  has  non-constant 
coefficients  whose  values  are  a  function  of  depth. 

When  the  eigenray  mode  is  selected,   rays  that  are  not 

refracted  or  reflected  to  within  yerror  of  the  receiver  are 

discarded. 
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III.   COMPUTER  SIMULATION  RESULTS 

A.  OVERVIEW 

The   computer   simulation   results  presented  in  this 
chapter  perform  phase  calculations  intended  to: 

-  show  that  ray  acoustics  can  be  used  to  evaluate 
the  phase  integral  producing  the  same  results  as 
direct  numerical  integration, 

-  show  that  the  two  methods  of  applying  ray  acoustics 
theory  (presented  in  Chapter  II)  produce  equivalent 
results  within  an  acceptable  tolerance, 

-  compare  speed  versus  accuracy  for  the  two  methods  for 
a  variety  of  test  cases,  and 

-  demonstrate  the  ability  to  identify  eigenrays  using 
either  of  the  ray  acoustics  methods. 

B.  DIRECT  INTEGRATION 

In  the  previous  chapter,  6  (y)  was  shown  as  the  solution 
to  the  phase  integral  in  the  depth  direction. 


8w(y)  =[Vk  (C)dC  (2.7) 

Jyo  y 


y 


The  straight-forward  evaluation  technique  is  the  brute 
force  approach  using  direct  numerical  integration.  The  phase 
integral  will  be  evaluated  using  direct  numerical  integration 
for  a  simple  test  case.  The  results  will  then  be  compared  to 
phase  values  obtained  using  ray  acoustics  calculating  phase 
values  from  travel  times. 
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Conditions  for  the  test  case  are: 

-  source  depth  of  10  meters, 

-  receiver  depth  of  100  meters, 

-  linear  sound-speed  profile  with  a  single,  constant 
gradient  of  0.016  sec   and 

-  four  values  of  radial  spatial  frequency  (FR  values) 
will  be  evaluated. 

For   each   FR  value,   tabular   results  will   show  the 

corresponding  ray  launch  angle  /3(y0).   The  relation  between 

the  FR  value  and  this  launch  angle  is  [Ref.  4]: 

)8(y0)  =  sin"1  [FR  *   c(y0)/f]  (3.1) 

where 

f  is  the  frequency  in  HZ,  and 

c(y0)  is  the  speed  of  sound  at  the  source  depth. 
Table  3 . 1  shows  the  results  produced  by  direct  numerical 
integration.   Tables  3.2  and  3.3  show  the  results  for  the 
Akima  cubic  spline/ordinary  differential  equation  solver 
(ODE  solver)  and  the  piecewise  linear  approximation  ray 
acoustics  techniques,  respectively. 

The  parameter  of  interest  is  THETAY  representing  the  phase 
change  in  the  depth  or  Y  direction.  Table  3.1  shows  these  6 
values  to  be  approximately  92,  86,  75,  and  56  radians  for  the 
four  FR  values  chosen.  Tables  3.2  and  3.3  contain  four  sets 
of  entries,  each  corresponding  to  an  FR  value.  These  entries 
begin  at  0  meters  range  and  end  at  the  range  corresponding  to 
100  meters  depth;  thus,  solving  the  phase  integral  for  the  6 
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TABLE   3.1 
PHASE    INTEGRAL   SOLUTION   BY  DIRECT   NUMERICAL    INTEGRATION 


INPUT  DATA  FOR  PHASE  INTEGRAL  EVALUATION 


F  •    25.0.0  MZ 
EABS  •  0.000000 
YREF  •     0.0  M 


VL  •    10.0  M 
EAEL  •  0.000100 
CYREF  ■  1500.0  H/SEC 


YU  •   100. 0  H 


RATIO  ■  o.m* 


C  •   0.016000  tl/SEC) 


EVALUATION  BY  1HSL10  ROUTINE  DODAG  WITH  IRULE 


FRICYCLES/M) 

BETAYO(DEG) 

RANGE    AT 

ANGLE    OF 

ARRIVAL 

TURNING    POINT 

TURNING    POINT 

THETAR(RAD) 

THETAYIRAD) 

l 

YU 

BETAY(DEG) 

DEPTHIH) 

RANGE (M) 

0.055:944859 

11.525 

18.56009 

11.536 

0.37555E'06 

0.45984E*06 

0.58408E'01 

0.«229167936E'02 

0.0665889718 

15.552 

59.25213 

23.576 

0.14090E'06 

0.21510E*06 

0.16425E'02 

0.8633840249E'02 

0.0998854576 

56.824 

67.45892 

36.866 

0.62692E'O5 

0.I2522E-06 

0.4I524E«02 

0.753785719IE*02 

0.  1551779455 

55.049 

119.60722 

53.122 

0.23574E«05 

0.70527E-05 

0. 10025E-05 

0.56S73684I5E'02 

TABLE    3.2 
PHASE    INTEGRAL    SOLUTION   BY   ODE    SOLVER 


ERAYS  *  F 

YO  ■    10.0  M 

FREOC/CHAX  «  0.16649  CYCLE3/M 

DLTFR  »   0.053294  CYCLES/H 

NOTEl  FRHIN  »  DLTFR 


PRTALL    =    T 

DEPTH    ■       100.0   M 

FREOC    •              250.0   HZ 

CHAX 

■       1501.6   H/SCC 

RATIO    •    0.9999 

FRHAX    ■    0.16647    CYCLES/H 

HFRS 

•      5 

RNG3TP    •       1 .0    M 

HRZRNG    •               150.0    H 

YR    • 

40.0    H 

FR  •   0.33294485881726E-01  CYCLE3/H 
BETAO  •      11.52457521170979  DEG 


RANGE(H)    DEPTH(M) 
0.0       10.00 
18.4      100.00 


TRVLTICEC)    BETAY(DEG)     BTEST        THETAT(RAD) 
0.000000      11.525    0.13318E-03    0. OOOOOOOE'OO 
0.061200      11.536    0.13318E-03    0.961 S2&3E-02 


THETAR(RAD) 
0. OOOOOOOE'OO 
0.3840846E'01 


THETAY(RAD) 
0. OOOOOOOE'OO 
0.9229168E»02 


FR  »   0.66S88971763452E-01  CYCLE3/H 
BETAO  .      23.55170283980132  DEG 


RANGE(M) 

DEPTHIH) 

TRVLT(SEC) 

0.0 

10.00 

0.000000 

39.3 

100.00 

0.065420 

BETAY(DEG)     BTEST 
23.552    0.26636E-03 
23.576    O.26636E-03 


THETAT(RAD) 
0. OOOOOOOE'OO 
0.1027611E«03 


THETAR(RAD)       THETAY(RAD) 
0. OOOOOOOE'OO     0. OOOOOOOE'OO 
0.1642273E'02     0.8633840E*02 


FR  «   0.99883457645178E-01  CYCLES/H 
BETAO  «      36.82440911644773  DEG 


RANGE (H) 

DEP1HIH) 

TRVLT(SEC) 

0.0 

10.00 

0.000000 

67.4 

100.00 

0.074932 

BETAY(DEG)     BTEST        THETAT(RAD) 
36.824    0.39953E-03    0 . O000O0OE '00 
36.866    0.39953E-03    0 . 1 1 7 7023E'03 


THETAR(RAD)        THETAY(RAD) 
0. OOOOOOOE'OO     0. OOOOOOOE'OO 
0.4232374E'02     0 . 753 7857E'02 


FR  =   0. 13317794552690E'00  CYCLES/H 
BETAO  «      53.04928575821315  DEG 


RANC-EIH)    DEPTH(M) 
0.0       10.00 
119.8       100.00 


TRVLT(SEC)    BETAY(DEG)     BTEST 
0.000000      53.049    0.53271E-03 
0.099859      53.122    0.55271E-05 


THETAT(RAD) 
0. OOOOOOOE'OO 
0.1568262E'03 


THETAR(RAD) 
0. OOOOOOOE'OO 
0. 1002S25E'03 


THETAYIRAD) 
0. OOOOOOOE'OO 
0.S657568E'02 
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TABLE    3.3 
PHASE    INTEGRAL    SOLUTION    BY    PIECEWISE    LINEAR   APPROXIMATION 


RAY  TRACING  USING  PIECEWISE  LINEAR  APPROXIMATION 
CASE  l*L<J 

NDATA  =   4  NUMBER  OF  GRADIENTS  =   3 

G(  1)  =   0.16000E-01  1/SEC  Gl  2)  *       0.16000E-01  1/SEC 
GC  J)  =   0.16000E-01  1/SEC 

ERAY3  =  F  PRTALL  =  T 

YO  =    10.0  M  DEPTH  =   100.0  M 
FREOC/CMAX  =  0.16649  CYCLE3/M     RATIO  »  0.9999 

DLTFR  =   0.033294  CYCLES/M  RNGSTP  *   1.0  M 
NOTE:  FRMIN  =  DLTFR 


FR  *   0.33294485881726E-01  CYCLES/M 
BETAO  •      11.52457521170980  DEG 

RANGE(M)    DEPTHCM)    TRVLTCSEC)    BETAY(DEG)     BTEST        THETAT(RAD)       THETARCRAD)       THETAY(RAD) 
0.0       10.00       0.000000      11.525    0.13318E-03    0 . 0000000E»00     0 . 00OO0OOE* 00     0 .  OOOOOOOE'OO 
18.  <J      100.00       0.061200      11.536    0.13318E-03    0 . 96 1 3253E»02     0 .  3840846E»0 1     0 .  9229 168E«02 


FREQC  «      250.0  HZ 

CMAX  »   1501.6  M/3EC 

FRMAX  ■  0.16647  CYCLES/M 

NFRS  ■   5 

HR2RNG  =•      150.0  M 

YR  =    40.0  M 

FR  :   0.66S88971763C5CE-01  CYCLE3/M 
EETAO  =      23.55170283980134  DEG 


RANGE(M) 

DEPTHCM) 

TRVLTCSEC) 

0.0 

10.00 

0.000000 

39.3 

100.00 

0.065420 

BETAY(DEG)     BTEST        THETAT(RAD)       THETARCRAD)       THETAY(RAD) 
23.552    0.26636E-03    0  .  OOOOOOOE'OO     0 . OOOOOOOE'OO     0 . 0000000E»00 
23.576    0.26636E-03    0 . 10276 1 1E*03     0 . 1 642273E*02     0 .863 3840E»02 


FR  -       0.99883457645178E-01  CYCLES/M 
BETAO  =      36.82440911644776  DEG 

RANGECM)    DEPTHCM)    TRVLTCSEC)    BETAYCDEG)     BTEST        THETATCRAD)       THETARCRAD)       THETAYCRAD) 

0.0       10.00       0.000000      36.824  0.39953E-03    0 . O000000E*00     0 . 0000000E* 00     0 .  O000000E*0O 

67.4      100.00       0.074932      36.866  0.39953E-03    0 .  1 1 77023E*03     0 .  4232374E»02     0  .  753 78S7E* 02 


FR  =   0. 13317794352690E*00  CYCLES/M 
EETAO  =      53.04928575821320  DEG 

RANGECM)    DEPTHCM)    TRVLTCSEC)    BETAYCDEG)     BTEST  THETATCRAD)       THETARCRAD)       THETAYCRAD) 

0.0       10.00       0.000000      53.049    0.53271E-03  0 . OOOOOOOE*00     0 . OOOOOOOE* 00     0 .  O0OOO0OE*O0 

119.8      100.00       0.099839      53.122    0.53271E-03  0 . 1568262E»03     0 . 100252SE' 03     0 .  5657 368E4 02 
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value  when  a  ray  travels  from  10  to  100  meters  depth.  The 
significance  of  these  three  tables  is  that  the  ray  acoustics 
methods  are  yielding  the  same  phase  values  0  as  the  direct 
integration  method. 

The  conclusion  to  be  derived  from  this  is  that  ray 
acoustics  theory  can  be  used  to  accurately  evaluate  the  phase 
integral . 

C.   VARYING  SOUND-SPEED  GRADIENTS 

Both  ray  acoustics  methods  will  now  be  applied  to  three 

different  test  cases.   A  zero  gradient,  positive  gradient 

-1  .  -1 

(+0.016  sec  ),  and  a  negative  gradient  (-0.016  sec  )  linear 

sound-speed  profile  will  be  used.  These  simulations  will  show 

that  both  methods  produce  equivalent  results  within  an 

acceptable  tolerance  for  a  variety  of  media. 

The  conditions  common  to  each  simulation  run  are: 

-  source  depth  of  10  meters, 

-  horizontal  range  of  1  kilometer, 

-  ocean  depth  of  100  meters, 

-  speed  of  sound  at  the  surface  of  1500  meters  per 
second,  and 

-  five   values   of   radial   spatial   frequency   (FR) 
evaluated. 

The  results  are  presented  in  Tables  3.4  through  3.9.  The 
phase  integral  solution  is  the  phase  change  in  the  Y  (depth) 
direction  listed  as  the  MODULO  TWOPI  THETAY  value  in  radians. 
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TABLE    3.4 
ZERO    GRADIENT,     PIECEWISE    LINEAR    SOLUTION 


RAY    TRACING    USING    PIECEWISE    LINEAR    APPROXIMATION 

CASE    1A 

NDATA  ■   4  NUMBER  OF  GRADIENTS  ■   3 

G(  1)  •   0.00000E-00  1/SEC    G(  2)  •   O.OOOOOE'OO  1/SEC 

Gl  S)  •   O.OOOOOE'OO  1/SEC 


VO  •    10.0  M 

FREOC/CMAX  ■  0.16667  CYCLES/M 
DLTFR  ■   O.OJ5330  CYCLES/M 
NOTE:  FRMIN  «  DLTFR 


DEPTH  • 

100.0  M 

RATIO  ■ 

0.9999 

RNGSTP  i 

•   1.0  M 

FREOC  ■  250.0  H2  CMAX  •   1500.0  M/SEC 

FPMAX  •  0.16645  CYCLES/M   NFRS  •   5 
HRZRNO  ■     1000.0  M 


R 

BETAO(DEG) 

DEPTH(M) 

TRVLT(SEC) 

BETAY(DEG) 

THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD) 

MODULO  TMOPI 
THETAY(RAD) 

1 

1 1 .536 

90.51 

3.333667 

168.464 

0.52365I1E»04 

0.2094186E*03 

0.5027093E-04 

0.5445951 

2 

23.576 

96.  44 

1.666833 

156.424 

0.2618256E»04 

0.4188371E*03 

0.2199419E«04 

0.3037135 

3 

36.866 

56.46 

1.111222 

143.134 

0.1745504E»04 

0.62825S7E»03 

0.1117248E»04 

S. 1243041 

< 

53.122 

39.79 

0.833417 

126.878 

0.1309128E*04 

0.8376743E«03 

0.4714536E«03 

0.2146886 

5 

89.190 

24.14 

0.666733 

89.190 

0.1047302E*04 

0.1047093E*04 

0.2094500E«00 

0.2094500 

TOTAL  CPU  TIME  ■  0  MIN  .  22.63  SEC 


TABLE    3.5 
ZERO    GRADIENT,     ODE    SOLVER    SOLUTION 


RAY  TRACING  USING  AKIMA  CUBIC  SPLINE  t  ODE  SOLVER 
CASE  2A 
NDATA  •   4 


YO  ■    10.0  M 

FREOC/CMAX  •  0.16667  CYCLES/H 
DLTFR  ■   0.033330  CYCLES/M 
NOTE:  FRMIN  ■  DLTFR 


DEPTH  •   100.0  M 

FREOC  •  250.0  HZ 

CMAX  • 

1500.0  M/SEC 

RATIO  •  0.9999 

FRMAX  •  0.16665  CYCLES/M 

NFRS  ■ 

5 

RNGSTP  ■   1.0  M 

HRZRNG  •     1000.0  M 

R 

BETAO(DEG) 

DEPTH(M) 

TRVLT(SEC) 

BETAY(DEG) 

THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD) 

MODULO  TWOPI 
THETAY(RAD) 

1 

11.536 

90.51 

3.333667 

168.464 

0.52!6511E»04 

0.2094186E«03 

0.5027093E*04 

0.5445951 

2 

23.576 

98.44 

1 .666833 

156.424 

0.26I82S6E«04 

0.4188371E«03 

0.2199419E*04 

0.3037135 

3 

36.866 

56.46 

1.111222 

143.134 

0.  1745504E«04 

0.6282557E«03 

0.1117248E'04 

5. 1243041 

4 

53.122 

39.79 

0.833417 

126.178 

0.  1309I28E«04 

0.8376743E«03 

0.4714536E*03 

0.2146886 

S 

89.190 

24.14 

0.666733 

8  9.190 

0.1047302E«04 

0.104709SE'04 

0.2094500E"00 

0.2094500 

TOTAL  CPU  TIME  »  15  MIN  .  52.38  SEC 
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TABLE    3.6 
POSITIVE    GRADIENT,     PIECEWISE    LINEAR    SOLUTION 


RAY  TRACING  USING  PIECEWISE  LINEAR  APPROXIMATION 

CASE  IB 

NDATA  •   4  NUMBER  OF  GRADIENTS  •   J 

G(  1  )  ■   0.16000E-OI  l/SEC    G(  2)  •   0.16000E-01  l/SEC 

G(  3)  ■   0.16000E-01  l/SEC 


VO  •    10.0  M 

FREOC/CMAX  •  0.1664"  CYCLES/M 
DLTFR  •   0.033294  CYCLES/M 
NOTE:  FRMIN  »  DLTFR 


DEPTH    •       100.0    M 

FREOC    •    250.0   H7 

CMAX    » 

1501 

RATIO    >    0.9999 

FRMAX    •    0.16647    CYCLES/M 

NFRS    • 

5 

RNGSTP    .       1.0    H 

HRZRNG    ■             1000.0    M 

R    1 

JETAO(DEG) 

DEPTH(M) 

TRVLT(SEC) 

BETAY(DEG) 

THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD) 

MODULO    TWOP1 
THETAY(RAD) 

1 

11.525 

87  .80 

3.335452 

168.466 

0.5236488E-04 

0.209I9S4E-03 

0.5027293E-04 

0.7442569 

2 

23.552 

96.99 

1  .666824 

156.425 

0.2«18241E*04 

0.4|83908E«03 

0.2I99850E-04 

0.7348999 

3 

.    56.824 

55.38 

1.111192 

143. 155 

0.I74S456E'04 

0.627S863E'0S 

0.  I1I7870E'04 

5.7462108 

4 

53.049 

38  .  73 

0.833378 

12* .927 

0. I309066E-04 

0.83678I7E>03 

0.4722847E'OS 

1.0457792 

5 

87. 3«3 

50.71 

0.667000 

87.974 

0. 1047721E-04 

O.I0459?7E«04 

0.1743798E«01 

1.7437978 

TOTAL  CPU  TIME  •  0  M1N  .  22.05  SEC 


TABLE    3.7 
POSITIVE    GRADIENT,    ODE    SOLVER   SOLUTION 


RAY  TRACING  USING  AKIMA  CUBIC  SPLINE  t  ODE  SOLVER 
CASE  2B 
NDATA  »   < 


YO  •    10.0  M 

FREOC/CMAX  •  0.16649  CYCLES/M 
DLTFR  ■   0.033294  CYCLES/M 
NOTE;  FRMIN  ■  DLTFR 


DEPTH    .       100.0    M 

FREOC    ■    250.0    HZ 

CMAX    • 

1501. t   M/SEC 

RATIO   •    0.9999 

FRMAX    ■    O.U447    CYCLES/M 

NFRS    ■ 

5 

RNGSTP    .       1 .0    M 

HRZRNG    ■             1000.0    M 

R 

BETAO(DEG) 

DEPTH(M) 

TRVLTISEC) 

BETAY(DEG) 

THETAT(RAD) 

THETARIRAD) 

THETAY(RAD) 

MODULO    TWOP1 
THETAY(RAD) 

1 

11.525 

87.80 

3.333652 

168.466 

0.5236488E*04 

0.209I954E-03 

0.5027293E-04 

0.7442693 

2 

23.552 

96.99 

1  .666824 

156.425 

0.26I824)E*04 

0.4183908E-03 

0.2I99850E'04 

0.7J49005 

S 

36.824 

55.38 

1  .111192 

143.155 

0.1745456E-04 

0.6275863E-03 

0.I1I7870E«04 

5.7462125 

4 

53.049 

38.73 

0.83:378 

126.927 

0.1309066E-04 

0.8S6781?E*03 

0.4722847E-03 

1.0457793 

5 

87.363 

50.71 

0.667000 

87.974 

0.  I04772IE-04 

0.1045977E«04 

0.1743798E-01 

1.7437978 

TOTAL  CPU  TIME  ■  16  MIN  .  24. 3C  SEC 
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TABLE    3.8 
NEGATIVE    GRADIENT,     PIECEWISE    LINEAR    SOLUTION 


RAY  TRACING  USING  PIECEWISE  LINEAR  APPROXIMATION 

CASE  1C 

NDATA  •   4  NUMBER  OF  GRADIENTS  •   5 

G(  1)  •  -0.16000E-01  I/SEC    G(  2)  •  -0.16000E-01  1/SEC 

G(    5)    '   -0.  H000E-01    l/SEC 


V0    •         10.0   M 

FREOC/CMAX  •  0.16667  CYCLES/M 
DLTFR  •   0.05SSJO  CYCLES/M 
NOTE:  FRMIN  •  DLTFR 


DEPTH  •   100.0  M 

FREOC  ■  250.0  HZ 

CMAX  ■ 

1500.0  M/SEC 

RATIO  •  0.9999 

FRMAX  «  0.16665  CYCLES/M 

NFRS  • 

5 

RNGSTP  •   1 .0  M 

HR2RNG  »     1000.0  M 

NFR  BETAO(DEG)  DEPTH(M)  TRVLT(SEC)  BETAY(DEG)   THETAT(RAD) 


11.555 
23.573 
36.861 
53.114 
88.835 


8  7.78 
96.98 
55.32 
38.63 
35.67 


I. 337240 
1  .668622 
1.112438 
0.834345 
0.667046 


168.475 
156.450 
143.160 
126.909 
88.224 


0.5242124E-04 
0.2621065E'04 
0.17474I4E-04 
0.  1310586E«04 
0.  1047793E»04 


THETAR(RAD) 

0.2094186E-03 
0.4188371E»03 
0.6282557E-03 
0.8376743E-03 
0.1047093E'04 


THETAY(RAD)  MODULO  TWOP  I 
THETAY(RAD) 
0.5032706E*04  6.1576516 
0.2202228E'04  5.1134653 
0.1119159E-04  0.7515713 
0.4729117E«03  1.6727921 
0.6999612E-00    0.6999612 


TOTAL  CPU  TIME  •  0  MIN  .  22.61  SEC 


TABLE    3.9 
NEGATIVE    GRADIENT,    ODE    SOLVER   SOLUTION 


RAY  TRACING  USING  AKIMA  CUBIC  SPLINE  t  ODE  SOLVER 
CASE  2C 
NDATA  *   4 


Y0  »    10. 0  M 

FREOC/CMAX  ■  0.16667  CYCLES/M 
DLTFR  •   0.033330  CYCLES/M 
NOTE:  FRMIN  •  DLTFR 


DEPTH  ■   100.0  M 

FREOC  •  250.0  HZ 

CMAX  » 

1500.0  M/SEC 

RATIO  •  0.9999 

FRMAX  '    0.16665  CYCLES/M 

NFRS  ■ 

5 

RNGSTP  •   1 .0  M 

HRZRNG  •     1000.0  M 

NFR  BETAO(DEG)  DEPTH(M)  TRVLT(SEC)  BETAY(DEG)   THETAT(RAD) 


11.535 
23.573 
36.861 
53.114 
88.835 


87.77 
96.98 
55.32 
38.6  3 
35.67 


3.337240 
1.668622 
1.112438 
0.834345 
0. 667046 


168.4  75 
156.450 
143.160 
126.909 
88.224 


0.5242124E«04 
0.2621065E-04 
0.1747414E«04 
0. I310586E»04 
0. 1047793E'04 


THETARIRAD1 

0.2094186E-03 
0.4188371E*03 
0.62825S7E-03 
0.8376743E'03 
0.1047093E«04 


THETAY(RAD) 

0.5032706E«04 
0.2202228E«04 
0.1119159E«04 
0.4729117E»03 
0.6999612E«00 


MODULO  TWOP  I 
THETAY(RAD) 
6.1576656 
3.1134683 
0.7515729 
1.6727921 
0.6999612 


TOTAL  CPU  TIME 
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Tables  3.4  and  3.5  show  the  results  for  the  zero  gradient 
case.  A  comparison  between  two  tables  shows  that  the  phase 
values  (MODULO  TWOPI  THETAY)  agree  perfectly.  Additionally, 
all  other  calculated  values  in  the  tables  agree.  The  only 
difference  is  the  TOTAL  CPU  TIME.  Both  methods  are  using  a 
conservative  range  step  (RNGSTP)  of  1  meter,  i.e.,  the 
incremental  increase  in  range  is  1  meter  in  the  propagation 
calculations.  For  these  conditions,  the  ODE  solver  is  taking 
over  3  0  times  longer  to  run  the  simulation. 

Tables  3.6  and  3.7  show  the  results  for  the  positive 
gradient  case.  A  comparison  between  these  two  tables,  as 
well  as  the  negative  gradient  results  shown  in  Tables  3.8  and 
3.9,  shows  acceptable  agreement  between  the  two  very 
different  calculation  methods.  The  only  discrepancies  in 
calculated  values  occurs  in  the  modulo  2n  0  phase  values 
needed  for  solving  the  phase  integral.  The  phase  calculation 
is  most  challenging  because  of  the  nature  of  phase  (being  a 
modulo  2n  function) .  The  discrepancies  for  both  the  positive 
and  negative  gradients  are  on  the  order  of  hundred-thousandths 
of  a  radian  or  less.  The  CPU  times  continue  to  follow  the 
pattern  seen  earlier  that  the  ODE  solver  requires  over  3  0 
times  longer  to  complete  a  simulation  run. 

The  following  ray  trace  plots,  Figures  3.1  through  3.6, 
correspond  to  Tables  3.4  through  3.9.  The  plots  assist 
interpretation  of  the  data,  but  as  expected  from  the  tabular 
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Figure  3.1   Ray  trace  corresponding  to  Table  3.4 
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Figure  3.2   Ray  trace  correseponding  to  Table  3.5 
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Figure  3.4   Ray  trace  corresponding  to  Table  3.7 
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Figure  3.5   Ray  trace  corresponding  to  Table  3.8 
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results,  no  variations  in  depth  or  angle  of  arrival  can  be 
detected. 

D.   CHALLENGING  TEST  CASES 

A  smooth  sound-speed  profile  in  the  form  of  a  half 
sinusoid  describes  the  medium  for  the  next  set  of 
simulations.  To  this  point,  the  agreement  in  accuracy  between 
the  two  methods  has  been  excellent;  however,  the  shape  of  a 
sine  curve  is  difficult  to  approximate  with  piecewise  linear 
segments.  While  this  will  be  a  challenging  test,  smoothly 
curved  profiles  do  occur  naturally.  The  following  set  of 
values  were  used: 

-  1500  m/sec  is  the  sound  speed  at  the  ocean  surface  and 
at  the  2000  meter  bottom,  and 

-  a  minimum  sound  speed  of  1475  m/sec  occurs  at  the  1000 
meter  depth. 

Three  sets  of  simulations  are  presented  using  5,  11,  and 
17  sound-speed  versus  depth  data  pairs.  The  simulation 
results  will  be  examined  to  quantify  the  number  of  ocean 
medium  sound-speed  samples  required  by  each  method  to  converge 
to  a  solution.  Also,  the  agreement  in  results  for  the  two 
methods  and  the  CPU  times  will  be  evaluated. 

As  the  number  of  equally  spaced  data  pairs  varies,  the 
medium  and  source  conditions  will  be: 

-  source  depth  of  1000  meters, 

-  horizontal  range  of  3  kilometers,  and 

five   values   of   radial   spatial   frequency   (FR) 
evaluated. 
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Tables  3.10,  3.11,  and  3.12  show  the  results  for  the  5, 
11,  and  17  data  pair  test  cases  for  the  piecewise  linear 
method.  In  comparing  the  first  two  of  these  tables,  large 
discrepancies  are  seen  in  the  final  depths  for  each  of  the 
five  radial  frequencies  evaluated — errors  as  large  as  seven 
meters.  Likewise,  travel  time  values  differ  by  hundredths  of 
a  second  meaning  that  modulo  2n  phase  values  show  no 
correlation  between  test  cases.  Unfortunately,  the  variations 
between  Tables  3.11  and  3.12  where  more  data  points  are  used, 
show  the  piecewise  linear  method  is  not  converging  to  a 
solution. 

Figures  3.7  through  3.9  correspond  to  Tables  3.10  through 
3.12.  While  numerical  phase  values  show  gross  disagreement, 
the  inability  to  detect  differences  in  the  graphical 
representations  illustrates  the  very  sensitive  nature  of  the 
phase  calculations. 

Tables  3.13  through  3.15  show  the  results  when  the  ODE 
solver  runs  the  simulations  for  5,  11,  and  17  data  pairs.  The 
5  data  pair  case  of  Table  3.13  again  shows  significant 
disagreement  with  the  Table  3.14  results  using  11  data  pairs. 
Final  depths  show  differences  of  as  much  as  a  meter  while 
travel  times  vary  by  milliseconds.  As  with  the  piecewise 
linear  method,  the  modulo  2w  phase  values  have  not  yet 
converged  to  a  solution.  Variations  of  up  to  2  radians  are 
seen. 
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TABLE    3.10 
FIVE    DATA    PAIRS,     PIECEWISE    LINEAR    SOLUTION 


RAY  TRACING  USING  PIECEWISE  LINEAR  APPROXIMATION 
CASE  1S05 


NDATA  >   5 

G(  1 )  «  -0.3S3S5E-01  1/SEC 

Gt  J)  ■   0.14645E-01  1/SEC 


NUMBER  OF  GRADIENTS  •   4 
G(  2)  •  -0. 14645E-01  1/SEC 
G(  4)  «   0.35355E-01  1/SEC 


ERAYS  ■  F 

YO  =  1000.0  M 

FREOC/CMAX  .  0.1*66?  CYCLES/M 

DLTFR  .   0.033530  CYCLES/M 

NOTE:  FRMIN  •  DLTFR 


PR  TALL  ■  F 

DEPTH  «  2000.0  M 

FREOC  ■     250.0  HZ 

CMAX  •   1500.0  M/SEC 

RATIO  •  0.9999 

FRMAX  «  0.16665  CYCLES/M 

NFRS  •   5 

RNG3TP  •   1.0  M 

HRZRNG  •     3000.0  M 

YR  ■    40.0  M 

R 

BETAO(DEG) 

DEPTH(M) 

TRVLT(SEC) 

BETAY(DEG) 

THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD) 

MODULO  TWOPI 
THETAY(RAD) 

1 

11.361 

144.74 

10.207186 

168.504 

0.1603341E-05 

0.6282557E«03 

0.1540515E*05 

5.0674943 

2 

23.160 

41.81 

5.103311 

156.449 

0.8016262E-04 

0.1256511E»04 

0.6759751E-04 

5.3264822 

3 

36.153 

1064.46 

3.402439 

36.180 

0.5344538E-04 

0.1884767E-04 

0.3459771E-04 

«'.0191220 

4 

51.868 

682.47 

2.554747 

127.902 

0.4012987E'04 

0.2513023E»04 

0. 1499964E»04 

4.5658515 

5 

79.494 

1510.35 

2.057830 

81.249 

0.3232433E-04 

0.3141278E-04 

0.9115404E*02 

S. 1894409 

TOTAL  CPU  TIME  ■  0  M1N  ,  24.95  SEC 


TABLE    3.11 
ELEVEN    DATA    PAIRS,     PIECEWISE    LINEAR    SOLUTION 


RAY    TRACING    USING    PIECEWISE    LINEAR    APPROXIMATION 
CASE    1S11 


NDATA 

•  11 

NUMBER  OF  GRADIENTS 

10 

G(  1) 

•  -0.38627E-01 

1/SEC 

G(  2)  •  -0. 34846E-01 

1/SEC 

G(  3) 

■  -0.27654E-01 

1/SEC 

G(  4)  .  -0.17755E-01 

1/SEC 

G(  5) 

■  -0.61179E-02 

1/SEC 

G(  6)  «   0.61179E-02 

1/SEC 

G(  7) 

■   0.17755E-01 

1/SEC 

G(  8)  =   0.27654E-01 

1/SEC 

G(  9) 

■   0.34846E-01 

1/SEC 

GUO)  «   0.38627E-01 

1/SEC 

ERAYS 

■  F 

PRTALL  •  F 

YO  «  1 

300.0  M 

DEPTH  •  2000.0  M 

FREOC  • 

250.0  HZ 

CMAX  • 

1500.0  M/SEC 

FREOC/ 

:HAX  »  0.16667 

CYCLES/M     RATIO  •  0.9999 

FRMAX  »  0 

16665  CYCLES/M 

NFRS  « 

5 

DLTFR 

■   0.033330  CYCLES/M 

RNG3TP  «   1  .0  M 

HRZRNG  ■ 

3000.0  M 

YR  ■ 

40.0  M 

NOTE:    FRMIN    »    DLTFR 


R 

SETAO(DEG) 

DEPTI 

HIM) 

TRVLT(SEC) 

JETA' 

riDEG 

)   THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD) 

MODULO  TWOPI 
THETAY(RAD) 

1 

11.341 

137 

50 

10.216764 

168 

506 

0. 1604846E«05 

0.6282557E«03 

0.1542020E-05 

1.2628446 

2 

23.160 

37 

97 

5.108073 

156 

449 

0.8023742E'04 

0.1256511E-04 

0.6767231E*04 

0.2400057 

5 

36.153 

1067 

39 

3.405626 

36 

165 

0.5349545E-04 

0.1884767E*04 

0.3464778E»04 

2.7430820 

4 

51.868 

679 

33 

2.557346 

127 

965 

0.4017070E*04 

0.2513023E-04 

0.1504048E«04 

2.3663691 

5 
OT 

79.494 

VL  CPU  TIME 

1520 
•  0 

40 
MIN 

2.060168 
.  25.36  SEC 

81 

354 

0.3236104E*04 

0.3141278E-04 

0.94825S0E*02 

0.5780165 
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TABLE  3.12 
SEVENTEEN  DATA  PAIRS,  PIECEWISE  LINEAR  SOLUTION 


RAY  TRACING  USING  PIECEWISE  LINEAR  APPROXIMATION 
CASE  IS17 


NDATA  i 

>  17 

NUMBER 

OF  GRADIENTS  > 

16 

G(  1)  > 

•  -o 

39018E- 

01 

1/SEC 

G(  2)  ■ 

-0 

37519E-01 

1/SEC 

G(  3)  < 

■  -o 

34577E- 

01 

1/SEC 

G(  4)  . 

-0 

30307E-01 

1/SEC 

G(  5)  ■■ 

>  -0 

24873E- 

01 

1/SEC 

G(  6)  ■ 

-0 

18482E-01 

1/SEC 

G(  7)  i 

.  -0 

11381E- 

01 

1/SEC 

G(  8)  ■ 

r  -o 

38429E-02 

1/SEC 

G(  9)  < 

>   0 

38429E- 

02 

1/SEC 

G(10)  > 

'   0 

11381E-01 

1/SEC 

G(  I  1  )  > 

■   0 

18482E- 

01 

1/SEC 

G(12)  < 

'   0 

24873E-01 

1/SEC 

G(I3)  < 

>   0 

30307E- 

■01 

1/SEC 

GC14)  < 

•   0 

34577E-01 

1/SEC 

G(15)  < 

>   0 

37S19E- 

■01 

1/SEC 

G(l«)  ■ 

>   0 

39018E-01 

l/SEC 

ERAYS  »  F 

Y0  «  "1000.0  M 

FREOC/CMAX  ■  0.16667  CYCLES/M 

DLTFR  •   0.033330  CYCLES/M 

NOTE:  FRMIN  •  DLTFR 


PRTALL  ■  F 
DEPTH  «  2000.0  M 
RATIO  '    0.9999 
RNGSTP  •   1.0  M 


FREOC  ■     250.0  HZ" 
FRMAX  •  0.16665  CYCLES/M 
HRZRNG  ■     3000.0  M 


CMAX  •   1500.0  M/SEC 

NFRS  •   5 

YR  »    40.0  M 


R 

BETAO(DEG) 

DEPTHCM) 

TRVLT(SEC) 

BETAY(DEG) 

THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD) 

MODULO  TWOPI 
THETAY(RAD) 

1 

11.341 

136.67 

10.217866 

168.506 

0.1605019E»05 

0.6282557E»03 

0.1542193E^05 

2.9937577 

2 

23.160 

37.53 

5.108620 

156.449 

0.8024602E*04 

0.I2S65I1E»04 

0.6768091E*04 

1.1002238 

3 

36.153 

1067.73 

3.405996 

36.160 

0.5350126E*04 

0.1884767E«04 

0.3465359E*O4 

3.3242727 

4 

51.868 

6  78.98 

2.557644 

127.973 

0.4017538E*04 

0.2513023E»04 

0.1504515E*04 

2.8338755 

5 

79.494 

1521.62 

2.060452 

81.322 

0.3236550E»04 

0.3141278E»04 

0.9527150E*02 

1.0237247 

TOTAL  CPU  TIME  »  0  MIN  .  25.94  SEC 
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Figure  3.7   Ray  trace  corresponding  to  Table  3.10 
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Figure  3.8   Ray  trace  corresponding  to  Table  3.11 
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Figure  3.9   Ray  trace  corresponding  to  Table  3.12 
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TABLE    3.13 
FIVE    DATA    PAIRS,     ODE    SOLVER 


RAY  TRACING  USING  AKIMA  CUBIC  SPLINE  t  ODE  SOLVER 
CASE  2S05 


NDATA  •  i 


ERAYS  •  F 

PRTALL  •  F 

YO  •  1000.0  M 

DEPTH  •  2000.0  H 

FREOC  ■      250.0  HZ 

FREOC/CMAX  •  t.Ult)  CYCLES/M 

RATIO  ■  0.9'99 

FRHAX  .  0.16665  CYCLES/H 

DLTFR  •   0.0JSS50  CYCLES/M 

RMGSTP  •   1 .0  H 

HRZRNG  •     5000. 0  M 

NOTE!  FRMIN  •  DLTFR 

CMAX  •   15  00.0  M/SEC 

NFRS  •   5 

YR  •    40.0  H 


R  BETAO(DEG) 

OEPTH(M) 

TRVLT(SEC) 

BETAYtDEG) 

THETATtRAD) 

THETAR(RAD) 

THETAY(RAD) 

MODULO  TWOPI 
THETAY(RAD) 

I    11.141 

154.85 

10.220257 

168.509 

0.  I605594E-05 

0.6282557E-0S 

O.I542S69E»05 

.  0.6666177 

2    25.160 

5«.5« 

S.10«825 

156.452 

0.8026491E>04 

0. 125651 IE-04 

0.  676998 OE«04 

2.9891864 

5    3  6.  155 

1068.49 

5.40(808 

56.157 

0.5551402E»04 

O.I884767E*04 

0.5466(55E*04 

4.5997015 

<    51.8(8 

678. 50 

2.558194 

127.975 

0.4018402E*04 

0  251502SE-04 

O.I50S579E-04 

!.  4978755 

5    79.494 

1521.82 

2.060498 

81.506 

0.525662SE*04 

0.3141278E*04 

0.9S544(8E<02 

1.0968988 

TABLE    3.14 
ELEVEN    DATA    PAIRS,    ODE    SOLVER 


RAY  TRACINO  USING  AKIMA  CUBIC  SPLINE  I  ODE  SOLVER 
CASE  2S11 


NDATA  •  II 


ERAYS  ■  F 

PRTALL  •  F 

Y0  •  1000.0  M 

DEPTH  •  2000.0  M 

FREOC  •      250.0  HZ 

FREOC/CHAH  .  0. 16667  CYCLES/M 

RATIO  ■  0.9999 

FRHAX  •  0.16665  CYCLES/M 

DLTFR  .   0.055550  CYCLES/M 

RNGSTP  ■   1.0  M 

HRZRNG  •     5000.0  M 

NOTEl  FRMIN  »  DLTFR 

CMAX  ■   1500. 0  M/SEC 

NFRS  •   5 

VR  •    40.0  M 


A  BETAO(DEG) 

DEPTH(M) 

TRVLT(SEC) 

BETAY(DEG) 

THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD) 

MODULO  THOP1 
THETAY(RAD) 

1    11.541 

156.09 

10.218625 

1(8.506 

O.I605158E*OS 

0.62825S7E»0S 

0.  I5425I2E-05 

4.1828982 

2    25.1(0 

57.25 

5.108997 

156.449 

0.8025I94E-04 

0  .  125651 1E*04 

0.4768685E*04 

1.4920(25 

5    56.155 

1067.97 

5.40(254 

56.157 

0.5550551E'04 

0.1884767E>04 

0.3465764E<04 

5.7288771 

4    51.868 

478.75 

2.557848 

127.578 

0.40I7858E'04 

0.251S025E-04 

0.1504855E«04 

1.1540595 

5    79.494 

1522.59 

2.0(0655 

81.512 

0.5256854E*04 

0.5141278E*04 

0.9555559E«02 

1.1078080 

TOTAL  CPU  TIME  •  48  MIN  .  14.55  SEC 


TABLE    3.15 
SEVENTEEN    DATA    PAIRS,     ODE    SOLVER 


RAY  TRACING  USING  AKIMA  CUBIC  SPLINE  I  ODE  SOLVER 
CASE  2517 


ERAYS  •  F 

PRTALL  •  F 

Y0  •  1000.0  M 

DEPTH  ■  2000.0  M 

FREOC  ■      250.0  HZ 

CMAX  •   1500.0  M/SEC 

FREOC/CMAX  .  0.16(67  CYCLES/M 

RATIO  ■  0.9999 

FRMAX  .  0.16665  CYCLES/M 

NFRS  ■   5 

DLTFR  .   0.055550  CYCLES/M 

RNGSTP  ■   1.0  M 

HRZRNG  •     5000.0  M 

YR  •    40.0  H 

NOTE:  FRMIN  ■  DLTFR 

R  BETAO(DEG) 

DEPTI- 

KM)  1 

RVLT(SEC)  1 

ETAY(DEO) 

THETAT(RAD) 

THETAR(RAD) 

THETAYIRAD) 

MODULO  TMOPI 
THETAYIRAD) 

1     11.541 

156 

15 

10.218580 

168.506 

0.1605151E«05 

0.6282557E*0S 

0.1542505E«05 

4.1156874 

2    25.140 

57 

25 

5. 108975 

156.449 

0.8025159E-04 

0.12S6511E-04 

0.47(B448E«04 

1.6570617 

5    54.155 

10(7 

96 

5.406259 

56.157 

0.S3S0S07E-O4 

0.1884767E-04 

0.5445740E«04 

3.7048568 

4    51.868 

6  78 

74 

2.557856 

12  7.978 

0.4017859E«04 

0.25I5025E-04 

0.  1504817E'04 

5.1555514 

5    79.494 

1522 

56 

2.0(0625 

81.517 

0.5256822E*04 

0.3I41278E»04 

0.9SS4574E>02 

1.2959589 

OTAL  CPU  TIME 

■  4 

MIN 

.  45.20  SEC 
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A  comparison  of  Tables  3.14  and  3.15  shows  that  phase 
calculation  results  converged  for  each  of  the  five  FR  values. 
The  increase  in  data  pairs  to  17  from  11  has  changed  the 
final  depth  values  by  only  a  few  hundredths  of  a  meter. 
Likewise,  the  very  sensitive  modulo  2T1  phase  values  show  a 
change  of  only  hundredths  of  a  radian  or  less.  The  ray  trace 
plots  for  Tables  3.13  through  3.15  are  provided  in  Figures 
3 . 10  through  3.12. 

The  CPU  times  for  these  three  tables  show  the  ODE  solver 
to  be  costly,  but  insensitive  to  increases  in  the  number  of 
data  pairs.  That  is,  approximately  the  same  amount  of  CPU 
time  is  required  regardless  of  the  number  of  data  pairs  used. 
The  CPU  time  for  the  piecewise  linear  method  also  proved 
fairly  insensitive  to  the  number  of  data  pairs  used  as  seen  in 
Tables  3.10  through  3.12. 

The  significant  findings  of  this  section  are: 

-  the  ODE  solver  can  perform  accurate  phase 
calculations  with  only  17  data  pairs  sampled  from  a 
2000  meter  deep  ocean  having  a  curved,  sinusoidal 
sound-speed  profile, 

-  the  piecewise  linear  approach  does  not  converge  to  a 
solution  using  the  17  data  pairs,  and 

-  the  differential  equation  solver  is  very  costly  to  use 

in  terms  of  CPU  time  versus  the  piecewise  linear  approach. 
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Figure  3.10   Ray  trace  corresponding  to  Table  3.13 
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Figure  3.11   Ray  trace  corresponding  to  Table  3.14 
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Figure  3.12   Ray  trace  corresponding  to  Table  3.15 
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E.   SPECIAL  INVESTIGATIONS 

Each  of  the  two  phase  calculation  methods  has  shown  one 
characteristic  weakness.  The  piecewise  linear  method  did  not 
converge  to  a  solution  in  the  sinusoidal  sound-speed  profile 
case  when  17  data  pairs  or  samples  were  provided.  In 
contrast,  the  ODE  solver  does  converge  to  a  solution,  but  is 
costly  in  terms  of  CPU  time.  This  section  will  attempt  to 
tailor  the  settings  for  each  method  to  compensate  for  its 
weaknesses. 

The  simulation  results  in  Tables  3.16  through  3.18  were 
produced  by  the  piecewise  linear  method  using  29,  55,  and  65 
data  pairs,  respectively.  All  medium  and  sound-speed  profile 
characteristics  remain  the  same  as  in  the  previous  section. 
As  the  number  of  data  pairs  increases,  the  6  values  are 
converging  to  the  Table  3.15  solutions  obtained  with  the  ODE 
solver.  Unfortunately,  the  phase  values  continue  to  show 
unacceptable  errors  of  tenths  and  hundredths  of  a  radian  for 
the  best  case  of  65  data  pairs. 

Table  3 . 19  shows  the  results  for  the  ODE  solver  when  the 
range  step  (RNGSTP)  parameter  is  increased  to  five  meters.  As 
explained  earlier,  the  range  step  is  the  differential 
equation  system  independent  variable.  All  calculations  prior 
to  Table  3 . 19  used  a  conservative  range  step  of  one  meter  for 
both  methods.  Table  3.19  compares  favorably  with  Table  3.15 
achieving  a  balance  between  CPU  time  and  accuracy.  While  the 
degradation  in  accuracy  is  only  microseconds  for  travel  time 
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TABLE  3.16 
TWENTY-NINE  DATA  PAIRS,  PIECEWISE  LINEAR  SOLUTION 


RAY  TRAC1N0  USING  PIECEWISE  LINEAR  APPROXIMATION 
CASE  1S29 


NDATA  > 

29 

NUMBER 

OF  GRADIENTS  ■ 

28 

C(  1)  • 

-0.39188E-01 

1/SEC 

G(  2)  • 

-0.3869SE-01 

1/SEC 

0(  J)  > 

-0.J7715E-01 

1/SEC 

0(  4)  • 

-0.S6262E-01 

1/SEC 

G(  S)  • 

-0.34352E-01 

1/SEC 

0(  «)  i 

-0.32010E-01 

1/SEC 

G(  7)  ■ 

-0.29266E-01 

1/SEC 

G(  8)  • 

■  -0.26154E-01 

1/SEC 

0(  9)  < 

-0.227I2E-01 

1/SEC 

G(10)  > 

>  -0.18986E-01 

1/SEC 

GUI)  • 

-0.IS020E-01 

1/SEC 

G(12)  • 

«  -0.10866E-01 

1/SEC 

6(13)  • 

r  -0.65745E-02 

1/SEC 

GU4)  • 

■  -0.22007E-02 

1/SEC 

0(15)  > 

■   0.22007E-02 

1/SEC 

0<1«)  < 

■   0.65745E-02 

1/SEC 

0(17)  ■ 

•   0.10866E-01 

1/SEC 

0(18) 

■   0.15020E-01 

1/SEC 

GU9)  < 

•   0.18986E-01 

1/SEC 

G(20) 

■   0.22712E-01 

1/SEC 

0(21) 

■   0.26154E-01 

1/SEC 

G(22) 

•   0.29266E-01 

1/SEC 

0(23) 

«   0.32010E-01 

1/SEC 

0(24) 

•   0.34352E-01 

1/SEC 

0(25) 

>   0.36262E-01 

1/SEC 

G(2i) 

■   0.37715E-01 

1/SEC 

0(27) 

■   0.38695E-01 

1/SEC 

G(28) 

■   0.39188E-01 

1/SEC 

ERAYS 

•  F 

PR 

TALL  ■  F 

YO  •  1 

)00.0  H 

DE 

»TH  ■  2000.0  H 

FREOC  ■ 

250.0  HZ 

CMAX  ■ 

1500.0 

M/SEC 

FREQC/ 

:HAX  >  0.16667 

CYCLES/M     RA 

TIO  •  0.9999 

FRMAX  •  0 

16665  CYC1.ES/M 

NFRS  ■ 

5 

DLTFR 

■   0.0SS330  CYCLES/M 

RM 

JSTP  •   1.0  M 

HR2RN0  ■ 

SOOO.O  M 

YR  ■ 

40.0  M 

NFR  BETAO(DEG)  DEPTH(M)  TRVLT(SEC)  BETAY(DEO)   THETAT(RAD) 


1      11 

341 

136 

31 

10 

218342 

1(8 

506 

0 

1(05093E*05 

2    23 

160 

37 

34 

5 

108857 

156 

449 

0 

8024973E»04 

3    3« 

153 

10(7 

.88 

3 

.406158 

36 

157 

0 

5350380E»04 

*    51 

8(8 

(78 

.82 

2 

557772 

127 

976 

0 

40177S9E*04 

5    79 

494 

1522 

.11 

2 

.0(0567 

81 

.320 

0 

.S236730E*04 

THETAR(RAD) 

0.6282557E*03 
0.12565UE*04 
0.1884767E*04 
0.2513023E»04 
0.3141278E*04 


THETAY(RAD)  MODULO  TNOPI 
THETAY(RAD 
0.1542268E»05  3.741S09 
0.6768462E*04  1.471084 
0.3465613E*04  3.577598 
0.1504716E»04  8.034874 
0.9545197E»02         1.204193 


TOTAL   CPU   TIME    ■    0   MIN    ,    26.46    SEC 
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TABLE  3.17 
FIFTY-FIVE  DATA  PAIRS,  PIECEWISE  LINEAR  SOLUTION 


RAY  TRACING  USING  PIECEWISE  LINEAR  APPROXIMATION 
CASE  1S55 


NDATA  ■ 

55 

NUMBER 

OF  GRADIENTS  ■ 

54 

G(  1)  ■ 

-0. 

S9248E-01 

1/SEC 

C(  2)  ■ 

-0. 

39115E-01 

1/SEC 

G(  3)  > 

-0. 

J88S0E-01 

1/SEC 

G(  4)  > 

-0. 

38453E-01 

1/SEC 

G(  S)  • 

-0. 

37926E-01 

1/SEC 

6(  t)  • 

-0. 

37271E-01 

1/SEC 

G<  7)  « 

-0. 

36490E-01 

1/SEC 

0(  8)  • 

-0. 

55584E-01 

1/SEC 

G(  9)  i 

-0. 

34541E-01 

1/SEC 

G(10)  • 

-0. 

33419E-01 

1/SEC 

GUI)  > 

-0. 

32K3E-01 

1/SEC 

GUI)  • 

-0. 

30800E-01 

1/SEC 

0(13)  ■ 

-0. 

29331E-01 

1/SEC 

G(14)  • 

-0. 

27764E-01 

1/SEC 

GU5)  > 

-0. 

2*103E-01 

1/SEC 

G(lt)  • 

-0. 

24353E-01 

1/SEC 

GU7)  > 

-0. 

22521E-01 

1/SEC 

6(18)  > 

-0 

20613E-01 

1/SEC 

Cll«)  • 

-0. 

18635E-01 

1/SEC 

G(20)  > 

-0. 

16594E-01 

1/SEC 

G(21)  > 

-0 

14497E-01 

1/SEC 

G(22)  i 

-0 

12350E-01 

1/SEC 

G(23)  • 

-0. 

1016ZE-01 

1/SEC 

0(24)  ■ 

-0 

79399E-02 

1/SEC 

0(25)  > 

-0 

56907E-02 

1/SEC 

G(24)  • 

>  -0 

34221E-02 

1/SEC 

0(27)  • 

-0 

11420E-02 

1/SEC 

G(28)  < 

>   0 

11420E-02 

1/SEC 

G(29)  • 

r   0 

34221E-02 

1/SEC 

0(30)  < 

•   0 

56907E-02 

1/SEC 

GUI)  ■ 

■   0 

79399E-02 

1/SEC 

0(32)  < 

>   0 

10162E-01 

1/SEC 

0(53)  < 

•   0 

12350E-01 

1/SEC 

0(34)  i 

•   0 

.14497E-01 

1/SEC 

0(35)  < 

>   0 

16594E-01 

1/SEC 

G(3«)  ■ 

•   0 

18«35E-01 

1/SEC 

0(57)  > 

.  6 

20*13E-01 

1/SEC 

0(38) 

•   0 

.22521E-01 

1/SEC 

0(39)  < 

•   0 

24353E-01 

1/SEC 

G(40) 

•   0 

.2610SE-01 

1/SEC 

6(41)  < 

•   0 

.27744E-01 

1/SEC 

0(42) 

•   0 

.293S1E-01 

1/SEC 

G(43) 

•   0 

.30800E-01 

1/SEC 

0(44) 

»   0 

.32163E-01 

1/SEC 

GK5) 

•   0 

.S3419E-01 

1/SEC 

G(4«) 

■   0 

.54561E-01 

1/SEC 

G(47) 

>   0 

.35586E-01 

1/SEC 

0(48) 

•   0 

.36490E-01 

1/SEC 

G(49) 

■   0 

.37271E-01 

1/SEC 

G(50) 

•   0 

.3792JE-01 

1/SEC 

6(51) 

•   0 

.384S3E-01 

1/SEC 

0(52) 

■   0 

.38850E-01 

1/SEC 

6(53) 

•   0 

.39115E-01 

1/SEC 

0(54) 

•   0 

.39248E-01 

1/SEC 

ERAYS  •  F 

YO  »  1000.0  M 

FREOC/CMAX    •    0.U447   CYCLES/M 

DLTFR    •      0.033330  CYCLES/M 

NOTE:    FRMIN   ■    DLTFR 


PRTALL    •   F 
DEPTH   •    2000.0   M 
RATIO   •    0.9999 
RNGSTP    «      1.0   M 


FREOC  ■     250.0  HZ 

CMAX  ■   1500.0  M/SEC 

FRMAX  -  0.1**65  CYCLES/M 

NFRS  •   5 

HRZRNG  ■     3000.0  M 

YR  •    40.0  M 

NFR   BETAO(DEG)    DEPTH(M)    TRVLT(SEC)    BETAY(DEG)      THETAT(RAD) 


11.341 
23.1*0 
3*. 153 

51.8(8 
79.494 


136.18 

37.28 

10*7.94 

«78.77 
1522.29 


10.218510 
5.108940 
3.40(215 
2.557817 
2.0*0*08 


1*8.50* 
15*. 449 

3*.  157 
127.977 

•1.317 


O.U05120E«05 
0.8025104E»04 
0.5350470E»04 
0.4017810E»04 
0.323*795E»04 


THETAR(RAD) 

0.*282557E«03 
0.125*511E»04 
0.18847*7E«04 
0.2513023E*04 
0.3141278E»04 


THETAY(RAD)  MODULO  TWOPI 
THETAY(RAD) 
0.1542294E*05  4.00*0532 
0.«7*8593E»04  1.6024783 
0.34«5702E»04  .  3.(*73828 
0.1504787E*04  3.1058334 
0.9551684E»02    1.2*90*45 


TOTAL  CPU  TIME  »  0  MIN  .  26.61  SEC 
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TABLE    3.18 
SIXTY-FIVE    DATA    PAIRS,    PIECEWISE    LINEAR   SOLUTION 


RAY  TRACING  USING  PIECEWISE  LINEAR  APPROXIMATION 
CASE  1S6S 


NDATA 
GI  1) 
Gf  3) 
G(  5) 
Gi  7) 
GI  9) 
GUI) 
G(13) 
GUS) 
G(17) 
G(l?) 
GCl  ) 
G(25) 
G(25) 
G(27) 
G<29) 
GCl) 
G(J3) 
G(35) 
G(37) 
G<39) 
S(«l) 
G<43) 
GC45) 
C<47) 
G(49) 
G(51) 
G(53) 
GC5S) 
G(57) 
G(59) 
GUI  ) 
G<63) 


65 

-0.39254E-01 

-0.38971E-01 

-0.38S12E-01 

-0.S7284E-01 

-0.35897E-01 

-0.34165E-01 

-0.32103E-01 

-0.29733E-01 

-0.27075E-01 

-0.24158E-01 

-0.21007E-01 

-0.17654E-01 

-0.14132E-01 

-0.10473E-01 

-0.67130E-02 

-0.28886E-02 

0.96364E-03 

0.48066E-02 

0.86032E-02 

0.12317E-01 

0.  15912E-01 

0.19354E-01 

0.22610E-01 

0.25647E-01 

0.28438E-01 

0.30955E-01 

0.33174E-01 

35073E-01 

36635E-01 

37844E-01 

38688E-01 

39160E-O1 


1/SEC 
I /SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
l/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 


NUMBER 

G(  2) 

G(  4  ) 

G(  «) 

G(  8) 

GI  10) 

GU2) 

GI  16) 

GI  It  ) 

GI  18) 

G(20) 

GI22) 

G(24) 

G(2t) 

G(28  ) 

G(30) 

G(32) 

G<34) 

G<36> 

GC8) 

G(40) 

G(42) 

GIU1 

GUI) 

G(48> 

G(50) 

G(52) 

GIS4) 

G(56) 

G(S8) 

G(tO) 

G<62) 

G(64) 


OF  GRADIENTS  ■ 
-0.39160E-01 
-0.38t88E-01 
-0.37844E-01 
-0.36635E-OI 
-0.35073E-01 
-0.33I74E-01 
-0.3095SE-0I 
-0.28438E-01 
-0.2S647E-01 
-O.22tl0E-0I 
-0.19354E-01 
-0.1S912E-01 
-0.12317E-01 
-0.8t032E-02 
-0.48066E-02 
-0.96364E-03 
0.28886E-02 
0.67130E-02 
10473E-01 
14132E-01 
17654E-01 
21007E-01 
24158E-01 
27075E-01 
29733E-01 
32103E-01 
34165E-01 
35897E-01 
37284E-01 
38312E-01 
38971E-01 


0.39254E-01 


64 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
t/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 
1/SEC 


ERAYS  =  F  PR  TALL  =  F 

Y0  =  1000.0  M  DEPTH  •  2000.0  M 

FREOC/CMAX  •  0 . 1«667  CYCLES/M  RATIO  ■  0.9999 

DLTFR  •   0.033330  CYCLES/M  RNGSTP  •   1.0  M 
NOTE:  FRMIN  •  DLTFR 


FREOC  ■  250.0  HZ  CMAX 
FRMAX  •  0.16665  CYCLES/M  NFRS 
HR2RNG  »    3000.0  M      YR  • 


1500.0  M/SEC 
5 
40.0  M 


R 

BETAO(DEG) 

DEPTH(M) 

TRVLTtSEC) 

BETAY(DEG) 

THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD) 

MODULO  TWOPI 
THETAY(RAD) 

1 

11.541 

136.17 

10.218528 

168.506 

0.1605123E»05 

0.6282557E*03 

0.1542297E*05 

«. 0340889 

2 

23.160 

37.27 

5.108949 

156.449 

0.8025118E»04 

0.12S6511E«04 

0.6768607E-04 

1.6163983 

3 

36.153 

1067.94 

3.406221 

36.157 

0.5350479E«04 

0.1884767E*04 

0.3465712E-04 

3.6769515 

4 

51 .868 

678. 76 

2.557822 

127.978 

0.4017818E'04 

0.2513023E»04 

0.150479SE-04 

3.1134251 

5 

79.494 

1522.31 

2  .060612 

81.317 

0.3236802E«04 

0.3141278E-04 

0.9552360E'02 

1.2758237 

TOTAL  CPU  TIME  »  0  MIN  ,  26.84  SEC 
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TABLE  3.19 
INCREASED  RANGE  STEP,  ODE  SOLVER  SOLUTION 


RAY  TRACING  USING  AKIMA  CUBIC  SPLINE  8  ODE  SOLVER 
CASE  2S17 


NDATA  ■  17 


ERAYS  •  F  PRTALL  ■  F 

YO  •  1000.0  M  OEPTH  .  2000.0  M     FREOC  •      250.0  HZ       CMAX  •   1500.0  M/SEC 

FREOC/CMAX  .  0.  1 66«7  CYCLES/M  RATIO  •  0.9999       FRMAX  .  0.16665  CYCLES/M   NFRS  •   5 

DLTFR  .   0.03J530  CYCLES/M  RNGSTP  »   5.0  M      HRZRNG  •     3000.0  M       YR  •    10.0  M 

NOTE:  FRMIN  .  DLTFR 


NFR  BETAO(DEG)  DEPTH(M)  TRVLT(SEC)  BETAY(DEG)   THETAT(RAD)     THETAR(RAD)  THETAY(RAD)   MODULO  TWOP I 

THETAY(RAD) 

1  11.311     136.09    10.218581   168.50*  0 . 1605131E«05  0 . 6282S57E«03  0. 1512306E-0S    <.  1219866 

2  23.160      37.25     5.108975   156.119  0  .8025159E-01  0 .  125651 1E«01  0. 6768618E'01    1.6570157 

3  36.153  1067.96  3.406239  36.157  0 .5350S07E«01  0 .  1881767E«01  0.  3165710E-01  3. 70X7770 
1  51.868  678.71  2.557836  127.978  0  .  10 l 7839E-01  0 . 25 l 3023E' 01  0 . 150181 7E«0«  3.1353500 
5    79.191    1522.36     2.060625    81.317  0 . 3236822E«01  0 . 3 11 12 78E «01  0 . 95513 79E«02    1. 2960153 

TOTAL  CPU  TIME  •  10  M1N  ,  3.53  SEC 


and  hundreds  or  thousandths  of  a  radian  for  modulo  2ir  8y,  the 
reduction  in  CPU  time  is  following  a  nearly  linear  relation 
to  the  increase  in  range  step,  that  is,  CPU  time  was  reduced 
approximately  by  one-fifth. 

F.   EIGENRAYS 

Tables  3.2  0  and  3.21  show  the  eigenrays  found  by  the  two 
ray  acoustics  methods.  Figures  3.13  and  3.14  correspond  to 
these  two  tables.  The  figures  show  the  sound-speed  profile  to 
be  that  of  a  SOFAR  channel.  This  channel  profile  was  chosen 
because  it  produces  interesting  ray  patterns,  and  it  commonly 
occurs  in  nature.  The  minimum  sound-speed  occurs  at  1000 
meters  depth  and  a  local  minima  occurs  at  the  ocean  surface. 
The  tables  show  that  each  method  found  the  same  eigenrays  at 
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TABLE    3.2  0 
EIGENRAYS    FOUND   USING    THE 
PIECEWISE    LINEAR   APPROXIMATION   METHOD 


RAY  TRACING  USING  PIECEWISE  LINEAR  APPROXIMATION 
CASE  ICH4 


NDATA  =   4 

Gt  1)  =   0.16000E-01  1/SEC 

G(  J)  =   0.17000E-01  1/3EC 


NUMBER  OF  GRADIENTS  «   5 
G(  2)  *  -0. I8889E-01  1/SEC 


ERAVS  =  T 

PR  TALL 

=  F 

YO  =   100. 0  M 

DEPTH  ■ 

2000.0  M 

FREOC  » 

250.0  HZ 

RNGSTP  »    5.0  M 

ANGSTP 

1.0  DEC 

1        HR2RNG 

•      200.0  H 

YR  =   250.0  H 

YERROR  •     5.0  M 

EIGENRAYS 

RAY  BETAO(DEG)  DEPIH(H) 

TRVLT(SEC) 

BETAY(DEG) 

THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD)   MO' 

I    53.000     251.11 

0. 16726? 

52.856 

0.2627462E«03 

0.I672491E-03 

0.95497I0E'02 

2    54.000     245.69 

0.165112 

53.856 

0.259356?E»03 

0.1694232E»03 

0.89?3371E'02 

S   150.000     246.73 

0.267010 

29.939 

0.419464?E»03 

0. IO47093E*O3 

0.3147556E«03 

MODULO  TWOPI 

THETAY(RAD) 

1 .2493187 

1 .9691119 

0.5963779 


TOTAL  CPU  TIME  *  0  MIN  ,  22. 47  SEC 


TABLE    3.21 

EIGENRAYS    FOUND   USING    THE 

ODE    SOLVER   METHOD 


RAY  TRACING  USING  AKIMA  CUBIC  SPLINE  8  ODE  SOLVER 
CASE  2CH4 


MDATA  *       4 

ERAYS  *    T 

PRTALL 

=  F 

YO  *   100.0  M 

DEPTH  : 

i  2000.0  M          FREOC  * 

250.0  H2 

RNGSTP  i    5.0  M 

ANGSTP 

I .0  DEG        HRZRNG 

•     200.0  H 

YR  =   250.0  M 

YERROR  =     5.0  M 

EIGENRAYS 

RAY  BETAO(DEG)  DEPIH(M) 

TRVLT(SEC) 

BETAY(DEG)   THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD)   MODULO  TWOPI 
THETAY(RAD) 

1    53.000     250.85 

0.167063 

52.930    0.2624216E«03 

0.1«72491E'03 

0.9517250E*02    0.924723? 

2    54.000     245.44 

0. 164912 

53.932    0.25<>0435E'03 

0. 16?4232E'03 

0.8°62032E'02    1.6557224 

3   150.000     246.56 

0.266343 

29.971    0.4191563E»03 

O.1047093E'03 

0.3144470E'03    0.2877732 

TOTAL  CPU  TIME  ■  24  MIN  ,  12. ?4  SEC 
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Figure  3.13   Ray  trace  plot  corresponding  to  Table  3.20 
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Figure  3.14   Ray  trace  plot  corresponding  to  Table  3.21 
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launch  angles  of  53,  54  and  150  degrees.  The  tables  highlight 
this  problem's  parameters  to  be: 

-  a  source  depth  y0  of  100  meters; 

-  a  receiver  depth  yr  of  2  50  meters; 

-  and  a  horizontal  range  (HRZRNG)  to  the  receiver  of  200 
meters. 

As  expected,  the  ODE  solver  runs  into  range  limitations 
imposed  by  CPU  time  requirements.  The  2  00  meter  example 
presented  consumed  over  24  minutes  of  CPU  time  using  the  ODE 
solver  versus  the  21  seconds  used  by  the  piecewise  linear 
method. 

Tables  3.22  and  3.23  (along  with  their  corresponding 
Figures  3.15  and  3.16)  show  solutions  to  long-range  eigenray 
search  problems.  These  two  examples  continue  to  use  a  source 
depth  y0  of  100  meters  and  a  receiver  depth  yr  of  250  meters. 
Table  3.22  shows  a  horizontal  range  (HRZRNG)  of  10  kilometers 
with  a  depth  error  yerror  at  the  receiver  of  2  meters.  The 
angle  step  size  (angstp)  between  each  ray  is  0.1°.  Table  3.22 
also  shows  that  three  eigenrays  were  found  with  launch  angles 
of  51.2°,  90.9°  and  129.5°.  The  CPU  time  for  this  10km  run  is 
shown  to  be  8.5  minutes. 
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TABLE    3.2  2 
LONG    RANGE    EIGENRAYS,     PIECEWISE    LINEAR    SOLUTION 


RAY  TRACING  USING  PIECEWISE  LINEAR  APPROXIMATION 
CASE  1CH4 


NDATA  =   4 

G(  1)  =   0.16000E-01  1/SEC 

G(  3)  =   0.17000E-01  1/SEC 


NUMBER  OF  GRADIENTS  =   3 
Gt  2)  =  -0.1888°E-01  1/SEC 


ERAYS  =  T 
YO  *       100.0  M 
RNGSTP  =    5.0  M 
YERROR  =     2.0  M 


PRTALL  =  F 
DEPTH  =  2000.0  M 
ANGSTP  =    0.1  DEG 


RAY  BETAO(DEG)  DEPTH(M)  TRVLT(SEC)  BETAY(DEG) 


1  51.200 

2  90.900 

3  129.500 


248.91 
250.77 
209.21 


8. 645958 
6.474996 

8.730378 


51.067 
86.356 

50.370 


FREOC  = 
HRZRNG 

EIGENRAYS 
THETAT(RAD) 

0.1358104E«05 
0.1048506E»05 
0.1371365E«05 


250.0  HZ 
10000.0  M 


THETAR(RAD) 

0.8160392E*04 
0.1046964E-05 
0.8079626E»04 


YR  =   250.0  M 


THETAY(RAD)   MODULO  TWOP1 

THETAY(RAD) 

0.5420648E«04    4.5420569 

0.1542213E'02    2.8557634 

0.5634020E«04        4.2855981 


TOTAL  CPU  TIME  =  8  MIN 


TABLE    3.2  3 
LONG   RANGE    EIGENRAYS,     PIECEWISE    LINEAR   SOLUTION 


RAY  TRACING  USING  PIECEWISE  LINEAR  APPROXIMATION 
CASE  1CH4 


NDATA  =   4 

C(  1)  =   0.16000E-01  1/SEC 

G(  3)  =   O.I7000E-01  1/SEC 


NUMBER  OF  GRADIENTS  =   3 
Gt  2)  =  -0. 18889E-01  1/SEC 


ERAYS  »  T 

PRTALL 

=  F 

V0  =   100.0  M 

DEPTH  = 

i  2000.0  M 

FREOC  = 

250.0  HZ 

RNGSTP  =    5.0  M 

ANGSTP 

1 .0  DEG        HRZRNG 

50000.0  M 

YR  =   250.0  M 

YERROR  =    15.0  M 

EIGENRAYS 

RAY  BETAO(DEG)  DEPTH(M) 

TRVLT(SEC) 

BETAY(DEG) 

THETAT(RAD) 

THETAR(RAD) 

THETAY(RAD) 

MODULO  TWOPI 

THETAY(RAD) 

1    40.000     236 

,69 

52.421823 

39.917 

0.8234401E»05 

0.3365291E»05 

0.4869109E-05 

2.6894315 

2   102.000     261 

,00 

34.440292 

102.535 

0.5409868E»05 

0.51210S7E»05 

0.2888117E*04 

4.1350528 

3   119.000     238 

,33 

38.520812 

60.820 

0.6050835E«05 

0.4579040E«05 

0.1<S7179SE»05 

2.7282732 

4   122.000     257 

.74 

39.734666 

122. 182 

0.6241507E«05 

0.4439925E*05 

0.1801581E«05 

1.9203099 

TOTAL  CPU  TIME  =  4  MIN  ,  27.61  SEC 
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Figure  3.15   Ray  trace  plot  corresponding  to  Table  3.22 
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Fl9Ure  3-16  *°>   —  ««=  corresponding  to 


Table  3.23 
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Table  3.23  shows  simulation  results  using  the  piecewise 
linear  method  with  a  horizontal  range  (HRZRNG)  of  50km. 
Parameters  that  affect  execution  time  such  as  angle  step  size, 
depths,  etc.,  match  the  simulation  presented  in  Table  3.21 
using  the  ODE  solver.  Table  3.21  shows  that  the  ODE  solver 
used  24  minutes  to  complete  a  200  meter  range  eigenvalue 
problem.  In  comparison,  Table  3.2  3  shows  that  the  piecewise 
linear  method  completed  the  50km  run  in  only  4   minutes. 

In  summary,  this  chapter  has  shown  that 

-  the  ODE  solver  performs  accurate  phase  calculations 
with  far  fewer  data  pair  samples  when  a  smooth 
sound-speed  profile  is  encountered; 

-  increasing  the  range  step  size  of  the  ODE  solver  to 
five  meters  reduces  the  CPU  time  cost  without 
signficantly  affecting  the  phase  solution  accuracies; 

-  the  piecewise  linear  method  can  run  long-range 
simulations  using  relatively  little  CPU  time;  and 

-  the  piecewise  linear  method  reguires  many  sound-speed 
data  pair  samples  for  an  accurate  solution  in  a  medium 
with  a  smooth,  sinusoidal  profile. 
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IV.   CONCLUSIONS  AND  RECOMMENDATIONS 

The  transfer  function  of  the  inhomogeneous  ocean  based  on 
the  WKB  approximation  requires  solving  a  phase  integral .  Ray- 
acoustics  theory  can  provide  solutions  to  this  phase  integral. 
The  simulations  performed  show  that  the  travel  times 
calculated  using  the  theories  of  ray  acoustics  can  be  used  to 
solve  the  phase  integral  avoiding  direct  numerical 
integration. 

Two  applications  of  ray  acoustics  produced  computer 
simulation  codes  which 

-  are  capable  of  solving  for  the  position,  travel  time 
and  phase  of  a  propagating  ray,  and 

-  have  very  different  advantages  and  costs. 

The  first  application  was  the  piecewise  linear  approximation. 
Sound-speed  versus  depth  data  pairs  sampled  from  the  ocean 
medium  were  connected  with  constant  gradient  linear  segments. 
Well-known,  closed  form  equations  form  the  mathematical  model 
for  sound  propagation.  The  solutions  are  low  cost  (in  terms 
of  CPU  time) ,  but  many  data  samples  are  required  for  accurate 
phase  solutions  for  arbitrary  sound-speed  profiles. 

In  contrast,  the  Akima  cubic  spline/ODE  solver  method  uses 
the  medium  samples  to  form  a  continuously  variable  sound-speed 
profile.  Accurate  phase  calculations  can  be  made  with  a 
minimum  of  sound-speed  versus  depth  data  pairs.  The 
disadvantage  in  using  the  ODE  solver  is  its  exorbitant  cost  in 
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terms  of  CPU  time.  Simulations  must  be  short  range  problems 
to  limit  computing  costs. 

The  ability  to  search  for  and  identify  eigenrays  was 
developed  for  each  of  the  two  ray  acoustics  theory 
applications.  Very  distinctive  characteristics  were  seen  for 
phase  calculations  in  terms  of  the  number  of  data  samples 
needed  to  assure  accurate  solutions.  The  identification  of 
eigenrays  is  an  easier  task  of  position  or  depth  computation. 
Here  the  distinctions  tend  to  disappear  with  the  two  methods 
arriving  at  the  same  solutions  for  eigenray  launch  angles. 
The  cost  or  CPU  time  required  continues  to  strongly  favor 
using  the  piecewise  linear  approach. 

These  findings  indicate  that  most  propagation  problems 
will  require  a  piecewise  linear  approach  for  computational 
efficiency.  Using  the  ODE  solver  would  severely  limit  the 
range  of  computer  simulations. 

Careful  sampling  of  the  ocean  media  is  required  to  obtain 
accurate  results  from  the  piecewise  linear  method.  Sufficient 
sound-speed  versus  depth  data  pairs  must  be  used  to  accurately 
represent  the  sound-speed  profile. 

This  thesis  developed  the  tools  to  quantify  the  strengths 

and  weaknesses  of  two  phase  computation  methods  in  a  variety 

of  media.    With  these   findings   in  mind,   future  work 

recommendations  are  to 

-  incorporate  each  phase  computation  technique  as  a 
module  in  the  larger  pulse  propagation  code;  and 
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-  run  pulse  propagation  simulations  to  compare  the 
received  pulse  shapes. 

If  this  comparison  shows  that  the  piecewise  linear 

approximation  produces  a  relatively  undistorted  received 

pulse,  it  is  an  efficient  solution  to  the  phase  computation 

problem. 


54 


REFERENCES 

1.  Ziomek,  Lawrence  J.,  Souza,  Luis  A.,  and  Campbell,  Peter 
R. ,  "Pulse  Propagation  in  a  Random  Ocean  —  A  Linear  Systems 
Theory  Approach,"  presented  at  OCEANS  '89,  Seattle,  WA,  pp. 
1211-1216,  September  1989. 

2.  Ziomek,  Lawrence  J.,  "Comments  on  the  Relationships 
Between  Linear  Systems  Theory  and  the  Free-space  Solution  of 
the  Inhomogeneous  Linear  Wave  Equation,"  J.  Acoust.  Soc.  Amer. 
88,  pp.  2027-2030,  October  1990. 

3.  Campbell,  Peter  R.  ,  An  Ocean  Medium  Pulse  Propagation 
Model  Based  on  Linear  Systems  Theory  and  the  WKB 
Approximation ,  Master's  Thesis,  Naval  Postgraduate  School, 
Monterey,  CA,  June  1989. 

4 .  Underwater  Acoustics  —  A  Linear  Systems  Theory  Approach , 
Academic  Press,  Orlando,  FL,  1985. 

5.  Lim,  C.K.,  Acoustic  Ray  Tracing,  Master's  Thesis,  Naval 
Postgraduate  School,  Monterey,  CA,  September  199  0. 


55 


INITIAL  DISTRIBUTION  LIST 


Defense  Technical  Information  Center 
Cameron  Station 
Alexandria,  VA  22304-6145 

Library,  Code  52 

Naval  Postgraduate  School 

Monterey,  California  93943-5002 

Chairman,  Code  EC 
Department  of  Electrical  and 
Computer  Engineering 
Naval  Postgraduate  School 
Monterey,  California  93943-5000 

Prof.  Lawrence  Ziomek,  Code  EC/Zm 
Department  of  Electrical  and 
Computer  Engineering 
Naval  Postgraduate  School 
Monterey,  California  93943-5000 

Prof.  Hung-Mou  Lee,  Code  EC/Lh 
Department  of  Electrical  and 
Computer  Engineering 
Naval  Postgradute  School 
Monterey,  California  93943-5000 

Dr.  Marshall  Orr 

ONR 

Code  11250A 

800  N.  Quincy  St. 

Arlington,  Virginia  22217 

CDR  Michael  L.  Whyms 
US  Naval  Academy 
Annapolis,  Maryland  23907 

Commanding  Officer 

ATTN:   LT  Daniel  L.  Devany 

Pearl  Harbor  Naval  Shipyard 

Box  400 

Pearl  Harbor,  Hawaii  96860-5350 


56 


Thesis 

D458715 

c.l 


Devany 

Sound  propagation  in 
the  inhomogeneous  ocean. 


Thesis 

D45&715  Devany 

c.l        Sound  propagation  in 
the  inhomogeneous  ocean. 


